function VarianceGammasTop23

addpath('basic_functions')
addpath('../Data')

Nboot = 1000;
Nfactors = 3;

[Return,Z,IsF,Nt,Chars,CharsNames,~,~,~] = package_data_all_greeks;
Return = Return.ret_daily;
[N,T,L] = size(Z);

f = [find(strcmp(CharsNames,'RV -- IV')) find(strcmp(CharsNames,'MarketCap')) find(strcmp(CharsNames,'Delta')) ...
     find(strcmp(CharsNames,'RV')) find(strcmp(CharsNames,'Max10')) ...
     find(strcmp(CharsNames,'Assets')) find(strcmp(CharsNames,'Stock price')) ...
     find(strcmp(CharsNames,'Moneyness')) find(strcmp(CharsNames,'Option price'))  find(strcmp(CharsNames,'Gamma')) ...
     find(strcmp(CharsNames,'Vega')) find(strcmp(CharsNames,'Turnover'))  find(strcmp(CharsNames,'Stock return')) ...     
     find(strcmp(CharsNames,'Rkurt')) find(strcmp(CharsNames,'RV -- MFvol'))  find(strcmp(CharsNames,'Debt')) ...     
     find(strcmp(CharsNames,'MFvol')) find(strcmp(CharsNames,'Leverage'))  find(strcmp(CharsNames,'IV ATM')) ...
     find(strcmp(CharsNames,'Stock return11')) find(strcmp(CharsNames,'IV slope'))  find(strcmp(CharsNames,'Profitability'))  find(strcmp(CharsNames,'Cash to asset'))];

Z = Z(:,:,[f 47]);
[N,T,L] = size(Z);

% construct W and X matrices that will be used for IPCA
W = NaN(L,L,T);
X = NaN(L,T);
for t = 1:T-1
    W(:,:,t) = squeeze(Z(IsF(:,t+1),t,:))'*squeeze(Z(IsF(:,t+1),t,:)) / Nt(t+1);
    X(:,t+1) = squeeze(Z(IsF(:,t+1),t,:))'*Return(IsF(:,t+1),t+1) / Nt(t+1);
    Return(~IsF(:,t+1),t+1) = NaN;
end
clear t;

% run IPCA
[G, F] = IPCA(X,W,Nt,Nfactors);
GammaBoot = NaN(Nboot,Nfactors,L);
for ell = 1:L
    [~,GammaBoot(:,:,ell)] = IPCA_char_pvalues(X,W,Nt,G,F,ell,Nboot);
end
save tables_outputs/GammaVariancesTop23 GammaBoot;

return